
;; restore,'./irdc_dist_model/bgps_rb3.sav',/ver
restore,'./irdc_dist_model/bgps_rb3_5647.sav',/ver
restore,'irdc_full_data.sav',/ver
s = read_bgps_csv(/ver)
n = n_elements(rb3)
print,'Number of RB3 sources: ',n


;; Make the TXT companion file
openw,lun,'./irdc_dist_model/examine_rb3.txt',/GET_LUN

printf,lun,'CNUM  GLON  GLAT  VLSR  CONTRAST  SIG_C  IMIR  SSUB  MORPH  IFACT'

FOR i=0L,n-1 DO BEGIN
   r = rb3[i]
   printf,lun,string(r.cnum,r.l,r.b,r.vlsr,r.c_meas,r.sig_c,r.i_mir,format=$
                     "(I4,'  ',F0.3,'  ',F0.3,'  ',F0.1,'  ',F0.3,'  ',"+$
                     "F0.3,'  ',F0.1,'    ')")
ENDFOR
close,lun
free_lun,lun


;; Make the figures

ellcol = 'Cyan'
ellth  = 5.
;; legcolor = 'Pale Green'
legcolor = 'WT3'
beamcol = 'Hot Pink'

last_glon = ''

FOR ll=0L, n-1 DO BEGIN
   
   cnum = string(rb3[ll].cnum,format="(I04)")
   message,'Plotting UBC for BGPS #'+cnum+'  ['+$
           string(ll+1,n,float(ll+1)/float(n)*100.,format=$
                  "(I4,' of ',I4,' (',F5.1,'%)]')"),/inf
   
   i = WHERE(s.cnum EQ rb3[ll].cnum,ni)
   
   ;; Read in the UBC frames
   ;; orig = readfits('local/ubc_glimpse_proc/nonsmooth'+cnum+'.fits',ohdr,/sil)
   smoo  = readfits('local/ubc_glimpse_proc/smooth30arc'+cnum+'.fits',shdr,/sil)
   imir  = readfits('local/ubc_glimpse_proc/Imir'+cnum+'.fits',ihdr,/sil)
   bgps  = readfits('local/ubc_glimpse_proc/bgps'+cnum+'.fits',bhdr,/sil)
   bsmoo = readfits('local/ubc_glimpse_proc/bsmoo'+cnum+'.fits',bshdr,/sil)
   label = readfits('local/ubc_glimpse_proc/label'+cnum+'.fits',lhdr,/sil)
   
   plot_wcs_axes,shdr,10,lonarr=xc,latarr=yc
   
   
   ;; Read in full GLIMPSE frame, if needed
   glon = string( round(s[i].glon_peak / 3.) * 3., format="(I03)")
   IF glon NE last_glon THEN BEGIN
      last_glon = glon
      undefine,gimg
      gimg = readfits('./local/glimpse_v3/starsub.'+glon+'00.I4.fits',ghdr)
      gsz  = size(gimg, /DIM)
   ENDIF
   
   
   gcirc,2,s[i].glon_peak,s[i].glat_peak,pf_irdc.l,pf_irdc.b,dis
   
   jj = WHERE(dis LE 240., njj)
   print,string(njj,format="('# of IRDCs: ',I0)")
   
   ;; Prep plotting
   myps,'./local/output/ubc_examine'+cnum+'.eps',xsize=11.5,ysize=10,/cmyk
   multiplot_xm, [2,2], mtitle='BGPS #'+cnum,ygap=0.03,xgap=0.35,/xmargin
   
   ;; Plot range
   pr = set_plot_range(alog(smoo))
   plotcs = 1.0
   
   ;;=============================================
   ;; Plot (a) -- STARSUB image
   
   ;; First do hextract to get smaller postage-stamp -- major speedup!
   adxy,ghdr,max(xc),min(yc),x0,y0
   adxy,ghdr,min(xc),max(yc),x1,y1
   hextract,gimg,ghdr,gps,gpsh,x0>0,x1<(gsz[0]-1),y0>0,y1<(gsz[1]-1),/silent
   
   plotmap,gps,gpsh,charsize=plotcs,ct=0,range=pr,/log
   
   FOR k=0L, njj-1 DO BEGIN
      j = jj[k]
      tvellipse,pf_irdc[j].min/3600.,pf_irdc[j].maj/3600.,$
                pf_irdc[j].l,pf_irdc[j].b,$
                -pf_irdc[j].pa_irdc,/data,thick=ellth,color=cgColor(ellcol)
   ENDFOR
   
   undefine,wtmask
   wtmask  = aperture_wtmask(gpsh, 40., s[i].glon_peak, s[i].glat_peak)
   plot_wcs_axes,gpsh,10,lonarr=xc,latarr=yc
   contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor(beamcol),/over
   al_legend,/top,/left,['(a)'],textcolor=legcolor,charsize=1.5,box=0
   
   multiplot

   ;;=============================================
   ;; Plot (b) -- BGPS Image w/ Label Contour
   
   undefine,bpr
   plot_wcs_axes,bhdr,10,lonarr=bxc,latarr=byc
   plotmap,bgps,bhdr,charsize=plotcs,ct=3,range=bpr,/log,min_dpi=100
   
   contour,label,bxc,byc,levels=[0.5],thick=5,color=cgColor('Green'),/over
   

   ;; NO IRDC ELLIPSES THIS FRAME
   ;; FOR k=0L, njj-1 DO BEGIN
   ;;    j = jj[k]
   ;;    tvellipse,pf_irdc[j].min/3600.,pf_irdc[j].maj/3600.,$
   ;;              pf_irdc[j].l,pf_irdc[j].b,$
   ;;              -pf_irdc[j].pa_irdc,/data,thick=ellth,color=cgColor(ellcol)
   ;; ENDFOR

   undefine,wtmask
   wtmask  = aperture_wtmask(shdr, 40., s[i].glon_peak, s[i].glat_peak)
   plot_wcs_axes,shdr,10,lonarr=xc,latarr=yc
   contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor(beamcol),/over
   al_legend,/top,/left,['(b)'],textcolor=legcolor,charsize=1.5,box=0
   
   multiplot
   
   ;;=============================================
   ;; Plot (c) -- Smoothed w/ BGPS flux contours
   plotmap,smoo,shdr,charsize=plotcs,ct=0,range=pr,/log,min_dpi=100
   
   ;; NO IRDC ELLIPSES THIS FRAME
   ;; FOR k=0L, njj-1 DO BEGIN
   ;;    j = jj[k]
   ;;    tvellipse,pf_irdc[j].min/3600.,pf_irdc[j].maj/3600.,$
   ;;              pf_irdc[j].l,pf_irdc[j].b,$
   ;;              -pf_irdc[j].pa_irdc,/data,thick=ellth,color=cgColor(ellcol)
   ;; ENDFOR
   
   undefine,wtmask
   wtmask  = aperture_wtmask(shdr, 40., s[i].glon_peak, s[i].glat_peak)
   plot_wcs_axes,shdr,10,lonarr=xc,latarr=yc
   contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor(beamcol),/over
   cgLoadct,40,/silent
   contour,bsmoo,bxc,byc,levels=[.025,.05,0.1,0.25,0.5],$
           c_colors=[1,2,3,4,5]*50,/overplot,thick=2
   cgLoadct,0,/silent
   
   al_legend,/top,/left,['(c)'],textcolor=legcolor,charsize=1.5,box=0
   
   multiplot
   
   ;;=============================================
   ;; Plot (d) -- I_MIR
   plotmap,imir,ihdr,charsize=plotcs,ct=0,range=pr,/log
   
   ;; NO IRDC ELLIPSES THIS FRAME
   ;; FOR k=0L, njj-1 DO BEGIN
   ;;     j = jj[k]
   ;;     tvellipse,pf_irdc[j].min/3600.,pf_irdc[j].maj/3600.,$
   ;;               pf_irdc[j].l,pf_irdc[j].b,$
   ;;               -pf_irdc[j].pa_irdc,/data,thick=ellth,color=cgColor(ellcol)
   ;;  ENDFOR
   
   undefine,wtmask
   wtmask  = aperture_wtmask(ihdr, 40., s[i].glon_peak, s[i].glat_peak)
   plot_wcs_axes,ihdr,10,lonarr=xc,latarr=yc
   contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor(beamcol),/over
   al_legend,/top,/left,['(d)'],textcolor=legcolor,charsize=1.5,box=0
   
   
   ;; Plot Colorbars
   cgLoadct,0,/silent
   cgcolorbar,range=exp(pr),/vertical,/right,position=[0.91,0.09,0.93,0.485],$
              /ylog,format="(F0.1)",title='GLIMPSE Intensity [MJy sr!u-1!n]',$
              charsize=plotcs

   cgLoadct,3,/silent
   cgcolorbar,range=exp(bpr),/vertical,/right,position=[0.91,0.545,0.93,0.94],$
              /ylog,format="(F0.3)",title='BGPS Flux Density [Jy beam!u-1!n]',$
              charsize=plotcs


   myps,/done,/mp
   
   ;; Clean up memory from loop
   undefine,gps,gpsh,wtmask,orig,smoo,imir

ENDFOR

END
